SUBROUTINE mayordecision(time,imun,iedu,ability,ipau,mwagein,iwealth,wealthin,ifu,fundsin,ipr,zzprin,iterm,pconsopt,qconsopt,zzpuopt,newsavopt,stealopt,utilopt,vemaxopt)
  USE commonvars
  USE Nelder_Mead
  USE simulated_anneal_mayor

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  REAL(8), INTENT(IN)  :: ability, mwagein, wealthin,fundsin,zzprin
  INTEGER, INTENT(IN)  :: time,imun,iedu,iwealth,ifu,ipr,iterm,ipau
  REAL(8), INTENT(OUT) :: pconsopt,qconsopt,zzpuopt,newsavopt,stealopt,utilopt,vemaxopt
  REAL(8)              :: a,b,fa,fb, avfundst, totresopt, emaxopt, netwealth, qcons_w_riv
  REAL(8)              :: stealopt0,newsavopt0,totresopt0,vemaxopt0,pconsopt0,qconsopt0,zzpuopt0,utilopt0,emaxopt0
  REAL(8)              :: stealopt1,newsavopt1,vemaxopt1,pconsopt1,qconsopt1,zzpuopt1,utilopt1,emaxopt1, Esteal, fine_fun
  REAL(8)              :: vemaxopt_sa, pconsopt_sa, qconsopt_sa, newsavopt_sa, stealopt_sa, utilopt_sa, emaxopt_sa, zzpuopt_sa
  INTEGER              :: iteration_sopt
  REAL(8), EXTERNAL    :: zeroin, focmayorTend

  INTEGER, PARAMETER   :: nn = 2, nn0 = 1
  REAL(8), ALLOCATABLE :: xx(:),xx0(:)

  ! Nelder-Mead             
  REAL(8)              :: object, simp, step(nn), stopcr, var(nn)
  INTEGER              :: iiiprint, iquad, maxf, nloop, ier
  LOGICAL              :: first
  ! Nelder-Mead 0             
  REAL(8)              :: object0, simp0, step0(nn0), stopcr0, var0(nn0)
  INTEGER              :: iiiprint0, iquad0, maxf0, nloop0, ier0
  LOGICAL              :: first0

  ! SIMULATED ANNEALING                                                                                                                       
  INTEGER, PARAMETER   :: neps1 = 1
  LOGICAL              :: ismax = .FALSE.

  EXTERNAL             :: objfcn_mayor, objfcn_mayor0, utilfun_mayor, prodfunc, dprodfunc

  iteration_sopt = 1

  IF (flag_simulation == 0) THEN
     fundst = per_funds*funds(ifu)
     zzprt = zzpr(ipr)
     netwealth = wealth(iwealth)
  ELSE
     fundst = fundsin
     zzprt = zzprin
     netwealth = wealthin
  END IF
!!$  netwealth = wealthin

  waget = mwagein
  imunt = imun
  iedut = iedu
  itermg = iterm
  ipaug = ipau
  timeg = time

  !High ability type
  ability_g = ability

  IF (time == Tend) THEN

     !Avaialble funds: restricting amount of available funds for stealing
     IF(imunt==3) THEN
        avfundst=omega3*fundst
     ELSE
        avfundst=omega1*fundst
     END IF
     !new savings must be zero in the last period
     newsavopt = wealth(w0)

     !total resources before stealing
     totrest = waget*lsmayor + R*netwealth - newsavopt

     !check if the the amount of available resources after stealing everything is positive
     !if not, there is no solution
     IF (totrest+avfundst <= subpcons) THEN

        pconsopt = subpcons + small_no
        qconsopt = subqcons + small_no
        zzpuopt = 0d0
        stealopt = 0d0 !fundst

        !We transform the level of public consumption to take into account the degree of rivalry
        !This is also the variable we output for public consumption
        !But, since in this case qconsopt = subqcons, we do not consider posizeg
!!$        qcons_w_riv = qconsopt/((popsizeg)**eta)
        qcons_w_riv = qconsopt

        CALL utilfun_mayor(pconsopt,qcons_w_riv,theta,utilopt)

        ! When we compute the value function we assign to the infeasibles a value that is a percentage of the lowest observed emax (in emaxfuntemp and emaxfun)
!!$        IF (flag_simulation == 0) utilopt = -infty

!!$        IF (flag_here == 647) write(*,*) '647 totres INF', totrest, waget*lsmayor, R*netwealth, newsavopt, finet, pconsopt,qconsopt,utilopt, subpcons

        RETURN

     END IF
     
     !bound used to constrain stealing from below
     boundt = MAX(-totrest+small_no,0d0)

     a = boundt
     b = avfundst

     fa = focmayorTend(a)
     fb = focmayorTend(b)

!!$     if (myrank == 0) write(*,'(A,4F30.10)') 'fa, fb', a, fa, b, fb

     IF(fa*fb < 0.0d0) THEN

!!$        write(*,'(A,4F12.4)') 'FA, FB, <0'

        !we look for the interior solution
        stealopt = zeroin(focmayorTend, a, b, 1d-10, 1d-10)
!!$        stealopt = boundt/(1.0d0+EXP(stealopt))+(fundst-small_no)*EXP(stealopt)/(1.0d0+EXP(stealopt))
     ELSE IF (fa < 0.0d0) THEN
        stealopt = a
     ELSE IF (fb > 0.0d0) THEN
        stealopt = avfundst - small_no
     END IF

     Esteal = paudit_g*stealopt + (1d0-paudit_g)*(stealopt-fine_fun(stealopt))

     IF (totrest + Esteal > 0d0) THEN
        pconsopt = totrest+Esteal
     ELSE
        pconsopt = small_no
     END IF
     
     !We use per-capita public inputs
     zzpuopt = (fundst - stealopt)/popsizeg
     IF (zzpuopt <= small_no) THEN
        qconsopt = subqcons
        GOTO 666
     END IF

     !The output qcons is per-capita; we account for the degree of rivalry in the utility function
     CALL prodfunc(imunt,iedut,zzpuopt,zzprt,qconsopt)

     !We are making sure that public consumption is not negative; In the simulations at the optimum it is never negative
     IF (qconsopt <= 0d0) qconsopt = subqcons

     !We transform per-capita public consumption in total public consumption
     !in emaxfuntemp and emaxfun we take into account the degree of rivalry by dividing by popsize**eta
666  qconsopt = qconsopt*popsizeg

     !We use the following variable only here to compute the utility level
     !We then take into account the degree of rivalry by dividing by popsizeg 
     qcons_w_riv = qconsopt/(popsizeg**eta)

     !The qcons is already in per-capita terms
     CALL utilfun_mayor(pconsopt,qcons_w_riv,theta,utilopt)

!!$     if (myrank == 0) write(*,'(a,8f20.10)') 'pcons, qcons_w_riv', pconsopt,qcons_w_riv,theta,utilopt,qconsopt,((popsizeg)**eta),popsizeg,eta

     IF (utilopt <  check_bignumber .OR. pconsopt <= 0.0d0) THEN
        write(*,*) 'util or pcons are negative in mayordecision 1', utilopt, pconsopt
        write(*,*) 'totrest,fundst,waget,lsmayor,netwealth,finet,stealopt', totrest,fundst,waget,lsmayor,netwealth,finet,stealopt
        write(*,*) 'a,b', a,b
        utilopt = bignumber
     END IF

!!$     IF (flag_here == 647) write(*,*) '647 totres', totrest, waget*lsmayor, R*netwealth, newsavopt, finet, pconsopt,qconsopt,utilopt

  ELSE

392  totrest = waget*lsmayor + R*netwealth

     !! Second, we calculate solution with stealing =0: solution0
     ! Define Initial Values steal = 0, saving = newsavopt - stealopt
     ALLOCATE(xx0(nn0))
     ! step sizes                   
     step0 = 1d0
     ! max # of function evaluations
     maxf0 = 10000
     ! print every iteration        
     iiiprint0 = -1
     ! stopping criterion           
     stopcr0 = 1d-8
     nloop0 = 2*nn0
     ! don't fit a quadratic surface
     iquad0 = 0
     ! accuracy 9 digits          
     simp0 = 1d-9
     ! declare this is first time 
     first0 = .true.

     vemaxoptg = -infty
     pconsoptg = subpcons
     qconsoptg = subqcons
     newsavoptg = netwealth
     stealoptg = 1d0
     utiloptg = -infty
     emaxoptg = -infty
     zzpuoptg = 0.0d0

     stealopt0 = 0.0d0
     totresopt = waget*lsmayor + R*netwealth + stealopt0
     IF (totresopt >= 0d0) THEN
        newsavopt = 0.1d0*totresopt
     ELSE
        newsavopt = 1.1d0*totresopt
     END IF
!!$     totresopt0 = waget*lsmayor + R*netwealth + stealopt
     xx0(1) = LOG((newsavopt - lbound_sav)/(totresopt-newsavopt))

     IF (isnan(xx0(1))) THEN
        WRITE(*,*) 'Sol 0 isnan',isnan(xx0(1)),xx0(1),newsavopt,lbound_sav,totresopt
        xx0(1) = LOG((netwealth - lbound_sav)/(totresopt-netwealth))
     END IF

     flag_sopt_sa = 0 
     CALL minim(xx0, step0, nn0, object0, maxf0, iiiprint0, stopcr0, nloop0, iquad0, simp0, var0, objfcn_mayor0, ier0)
     
     stealopt0 = 0d0
     totresopt0 = waget*lsmayor + R*netwealth + stealopt0
     newsavopt0 =  (totresopt0*EXP(xx0(1)) + lbound_sav)/(1.0d0 + EXP(xx0(1)))
     
     IF (ier0 > 0) THEN
        write(*,'(a,1i4,3f20.10,1f50.10)') 'minim0 not a solution', ier0, xx0, stealopt0,newsavopt0,object0
        write(*,'(a,8I3,1f8.4)') 'state no sol:t,mun,ed,ipau,au,iw,ifu,pr,term,ability',time,imun,iedu,ipau,iwealth,ifu,ipr,iterm,ability
     END IF

     vemaxopt0 = vemaxoptg  
     pconsopt0 = pconsoptg  
     utilopt0 = utiloptg  
     emaxopt0 = emaxoptg  

     !We use per-capita public inputs
     zzpuopt0 = (fundst - stealopt0)/popsizeg
     IF (zzpuopt0 <= small_no) THEN
        qconsopt0 = subqcons
        GOTO 667
     END IF

     !The output qcons is per-capita; we account for the degree of rivalry in the utility function
     CALL prodfunc(imun,iedu,zzpuopt0,zzprt,qconsopt0)

     !We transform per-capita public consumption in total public consumption
     !in emaxfuntemp and emaxfun we take into account the degree of rivalry by dividing by popsize**eta
667  qconsopt0 = qconsopt0*popsizeg

!!$     IF(flag_here==642) WRITE(*,'(a,1F20.12,2I3)') 'res,iterm,iwealth',totrest,iterm,iwealth
!!$     IF(flag_here==651) WRITE(*,'(a,7F20.14)') 'SolMinim0',pconsopt0,qconsopt0,newsavopt0,stealopt0,utilopt0,emaxopt0,vemaxopt0
     DEALLOCATE(xx0)

     !xx is the initial condition on input and the solution on output
     ALLOCATE(xx(nn))
!!$101  ALLOCATE(xx(nn))

     IF (iteration_sopt == 1 .OR. (iteration_sopt == 2 .AND. flag_nan_minim == 1)) THEN
!!$     IF (iteration_sopt == 2 .AND. flag_nan_minim == 1) THEN

        iteration_sopt = 1
        flag_nan_minim = 0

        ! step sizes                   
        step = 1d0
        ! max # of function evaluations
        maxf = 10000
        ! print every iteration        
        iiiprint = -1
        ! stopping criterion           
        stopcr = 1d-8
        nloop = 2*nn
        ! don't fit a quadratic surface
        iquad = 0
        ! accuracy 9 digits          
        simp = 1d-9
        ! declare this is first time 
        first = .true.

        !totrest = waget*lsmayor + R*netwealth
        totresopt = waget*lsmayor + R*netwealth + stealopt0
!!$        xx(1) = LOG((stealopt0+small_no)/(fundst-(stealopt0+small_no)))
        xx(1) = LOG((stealopt0+small_no+0.05d0*fundst)/(fundst-(stealopt0+small_no)))

        if (ABS(totresopt-newsavopt0) < 0.0001d0) WRITE(*,'(a,4F20.10)') 'divided by zero', totresopt-newsavopt0, totresopt, newsavopt0, lbound_sav
        if ((newsavopt0 - lbound_sav)/(totresopt-newsavopt0) < 0.0001d0) THEN
           WRITE(*,'(a,7F30.15)') 'close to zero', (newsavopt0 - lbound_sav)/(totresopt-newsavopt0), totresopt, newsavopt0, lbound_sav, waget*lsmayor, netwealth, stealopt0
!!$           xx(5) = 1d0
        END if
        xx(2) = LOG((newsavopt0 - lbound_sav)/(totresopt-newsavopt0))

        ismax = .FALSE.
        vemaxoptg = -infty
        pconsoptg = subpcons
        qconsoptg = subqcons
        newsavoptg = netwealth
        stealoptg = 1d0
        utiloptg = -infty
        emaxoptg = -infty

!!$        IF(flag_here==246) WRITE(*,'(a,9F20.10)') 'mayordec b 1', xx(1), xx(2), stealopt0, newsavopt0, totresopt, fundst, pconsopt0, utilopt0, emaxopt0

        flag_sopt_sa = 0 
        CALL minim(xx, step, nn, object, maxf, iiiprint, stopcr, nloop, iquad, simp, var, objfcn_mayor, ier)

!!$        IF(flag_here==246) WRITE(*,'(a,4F20.10)') 'mayordec a 1', xx(1), xx(2), (EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst, ((totresopt+(EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst)*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2)))

        stealopt_sa = stealoptg
        newsavopt_sa = newsavoptg
        vemaxopt_sa = vemaxoptg
        pconsopt_sa = pconsoptg
        utilopt_sa = utiloptg
        emaxopt_sa = emaxoptg

!!$        ! WE find initial values with SA
!!$        ALLOCATE(ubp_sopt(nn),lbp_sopt(nn))
!!$
!!$        lbp_sopt = (/  0.0d0, -500000000d0/norm /)
!!$        ubp_sopt = (/ fundst, totrest + 0.85d0*fundst /)
!!$        vm = (/ fundst, 1000000d0/norm /)
!!$     
!!$        xx(1) = stealopt0
!!$        xx(2) = newsavopt0
!!$
!!$        flag_sopt_sa = 1
!!$        CALL sa(nn, xx, ismax, rt, eps, ns, nt, neps1, maxevl, lbp_sopt, ubp_sopt, c, iiprint, iseed1, iseed2, t, vm, xopt, fopt, nacc, nfcnev, nobds, ier)

!!$        !Translating solution into local vars
!!$        vemaxopt_sa = vemaxoptg  
!!$        pconsopt_sa = pconsoptg  
!!$        newsavopt_sa = xx(2)
!!$        stealopt_sa = xx(1)
!!$        utilopt_sa = utiloptg  
!!$        emaxopt_sa = emaxoptg  



        !We use per-capita public inputs
        zzpuopt_sa = (fundst - stealopt_sa)/popsizeg
        IF (zzpuopt_sa <= small_no) THEN
           qconsopt_sa = subqcons
           GOTO 668
        END IF

        !The output qcons is per-capita; we account for the degree of rivalry in the utility function
        CALL prodfunc(imun,iedu,zzpuopt_sa,zzprt,qconsopt_sa)

        !We transform per-capita public consumption in total public consumption
        !in emaxfuntemp and emaxfun we take into account the degree of rivalry by dividing by popsize**eta
668     qconsopt_sa = qconsopt_sa*popsizeg

!!$        DEALLOCATE(ubp_sopt,lbp_sopt)

     END IF

     !! First we calculate solution with no constraints on  stealing and savings: solution1

     !Now use Minim with SA solution as starting values
     ! step sizes                   
     step = 1d0
     ! max # of function evaluations
     maxf = 10000
     ! print every iteration        
     iiiprint = -1
     ! stopping criterion           
     stopcr = 1d-8
     nloop = 2*nn
     ! don't fit a quadratic surface
     iquad = 0
     ! accuracy 9 digits          
     simp = 1d-9
     ! declare this is first time 
     first = .true.

     !totrest = waget*lsmayor + R*netwealth
23   totresopt = waget*lsmayor + R*netwealth + stealopt_sa
     IF (fundst-stealopt_sa > small_no) THEN
        xx(1) = LOG(stealopt_sa/(fundst-stealopt_sa))
     ELSE
        xx(1) = LOG((stealopt_sa*0.5d0)/(fundst-stealopt_sa*0.5d0))
     END IF
     IF (stealopt_sa < small_no) THEN
        xx(1) = LOG((stealopt_sa+small_no)/(fundst-(stealopt_sa+small_no)))
     END IF
     xx(2) = LOG((newsavopt_sa - lbound_sav)/(totresopt-newsavopt_sa))

!!$23   steal_temp = 0.1*fundst
!!$     totresopt = waget*lsmayor + R*netwealth + steal_temp
!!$     xx(1) = LOG((stealopt0)/(fundst-steal_temp))
!!$     IF (steal_temp < small_no) THEN
!!$        xx(1) = LOG((steal_temp+small_no)/(fundst-(steal_temp+small_no)))
!!$     END IF
!!$     xx(2) = LOG((newsavopt0 - lbound_sav)/(totresopt-newsavopt0))

!!$     if (flag_here == 890) WRITE(*,'(a,10f25.10)') 'xx', xx, totresopt, waget*lsmayor, netwealth, newsavopt0, stealopt0, lbound_sav, newsavopt - lbound_sav, totresopt-newsavopt0

     IF (isnan(xx(1)) .OR. isnan(xx(2))) THEN
        iteration_sopt = 2
        flag_nan_minim = 1
        WRITE(*,*) 'Sol isnan BIG PROBLEM',isnan(xx(1)),isnan(xx(2))
        WRITE(*,'(a,13f25.10)') 'BIG PRoblem xx', xx, totresopt, newsavopt_sa, stealopt_sa, lbound_sav, newsavopt_sa - lbound_sav, totresopt-newsavopt_sa, fundst, fundst-stealopt_sa, pconsopt_sa, utilopt_sa, emaxopt_sa
!!$        deallocate(xx)
!!$        GO TO 101
     END IF

     vemaxoptg = -infty
     pconsoptg = subpcons
     qconsoptg = subqcons
     newsavoptg = netwealth
     stealoptg = 1d0
     utiloptg = -infty
     emaxoptg = -infty
     zzpuoptg = 0.0d0

!!$     IF(flag_here==246) WRITE(*,'(a,9F20.10)') 'mayordec b 2', xx(1), xx(2), stealopt_sa, newsavopt_sa, totresopt, fundst, pconsopt_sa, utilopt_sa, emaxopt_sa

     flag_sopt_sa = 0 
     CALL minim(xx, step, nn, object, maxf, iiiprint, stopcr, nloop, iquad, simp, var, objfcn_mayor, ier)

!!$     IF(flag_here==246) WRITE(*,'(a,4F20.10)') 'mayordec a 2', xx(1), xx(2), (EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst, ((totresopt+(EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst)*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2)))

!!$     IF(flag_here==178) WRITE(*,'(a,1i3,7F24.14)') 'Minim sol', ier, stealoptg, (EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst, newsavoptg, ((waget*lsmayor + R*netwealth + stealopt1)*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2))), vemaxoptg, utiloptg, emaxoptg

     IF (vemaxoptg > vemaxopt_sa .AND. ier <= 0) THEN
        stealopt1 = stealoptg !(EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst
!!$        totresopt1 = waget*lsmayor + R*netwealth + stealopt1
        newsavopt1 = newsavoptg !(totresopt1*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2)))
        vemaxopt1 = vemaxoptg  
        pconsopt1 = pconsoptg  
        utilopt1 = utiloptg  
        emaxopt1 = emaxoptg  

        !We use per-capita public inputs
        zzpuopt1 = zzpuoptg
        IF (zzpuopt1 <= small_no) THEN
           qconsopt1 = subqcons
           GOTO 669
        END IF

        !The output qcons is per-capita; we account for the degree of rivalry in the utility function
        CALL prodfunc(imun,iedu,zzpuopt1,zzprt,qconsopt1)

        !We transform per-capita public consumption in total public consumption
        !in emaxfuntemp and emaxfun we take into account the degree of rivalry by dividing by popsize**eta
669     qconsopt1 = qconsopt1*popsizeg

     ELSE
        stealopt1 = stealopt_sa
        newsavopt1 = newsavopt_sa
        vemaxopt1 = vemaxopt_sa
        pconsopt1 = pconsopt_sa
        utilopt1 = utilopt_sa
        emaxopt1 = emaxopt_sa
        zzpuopt1 = zzpuopt_sa
        qconsopt1 = qconsopt_sa
     END IF

!!$     IF(flag_here == 246) WRITE(*,'(a,9F20.10)') 'sol sa', vemaxopt_sa, stealopt_sa, newsavopt_sa, pconsopt_sa, qconsopt_sa, zzpuopt_sa, utilopt_sa, emaxopt_sa, paudit_g
!!$     IF(flag_here == 246) WRITE(*,'(a,9F20.10)') 'sol mi', vemaxoptg, stealoptg, newsavoptg, pconsoptg, qconsoptg, zzpuoptg, utiloptg, emaxoptg, paudit_g

!!$     IF (ier > 0) THEN
!!$
!!$        write(*,'(a,1I3,7f20.10)') 'minim not a solution mayor', ier, xx, (EXP(xx(1))/(1d0 + EXP(xx(1))))*fundst,(totresopt1*EXP(xx(2)) + lbound_sav)/(1.0d0 + EXP(xx(2))), object, vemaxoptg, vemaxopt_sa


!!$        WRITE(*,'(a,2f30.10,7F20.10,2f30.10)') 'minim not a solution mayor',totresopt,netwealth,pconsoptg,qconsoptg,qconsoptg/((popsizeg)**eta),fundst,newsavoptg,stealoptg,utiloptg,emaxoptg,vemaxoptg
!!$        WRITE(*,'(a,2f30.10,7F20.10,2f30.10)') 'sa          solution mayor',totresopt,netwealth,pconsopt_sa,qconsopt_sa,qconsopt_sa/((popsizeg)**eta),fundst,newsavopt_sa,stealopt_sa,utilopt_sa,emaxopt_sa,vemaxopt_sa
!!$        write(*,'(a,10I3,1f8.4)') 'state no sol:t,mun,ed,ipau,au,iw,ifu,pr,term,ability',time,imun,iedu,ipau,iwealth,ifu,ipr,iterm,ability
!!$        iteration_sopt = 2
!!$        flag_nan_minim = 1
!!$        flag_here = 691
!!$        deallocate(xx)
!!$        GO TO 101

!!$     END IF

!!$     IF (vemaxoptg > 100000d0) THEN
!!$        WRITE(*,'(a,2f30.10,7F20.14)') 'vemaxoptg > 0',totresopt,netwealth,pconsoptg,qconsoptg,newsavoptg,stealoptg,utiloptg,emaxoptg,vemaxoptg
!!$        flag_here = 691
!!$        GO TO 23
!!$     END IF
     
     DEALLOCATE(xx)

!!$     IF(flag_here==429) WRITE(*,'(a,7F20.14)') 'SolMinim1',pconsopt1,qconsopt1,newsavopt1,stealopt1,utilopt1,emaxopt1,vemaxopt1
!!$     IF(flag_here==429) WRITE(*,'(a,7F20.14)') 'SolMinim0',pconsopt0,qconsopt0,newsavopt0,stealopt0,utilopt0,emaxopt0,vemaxopt0
!!$     IF(flag_here==429) WRITE(*,*) ''

     !! Choosing between the two maximum
     IF(vemaxopt0 > vemaxopt1 .OR. stealopt1 < 0.01d0) THEN
        stealopt = stealopt0
        newsavopt = newsavopt0
        vemaxopt = vemaxopt0  
        pconsopt = pconsopt0
        !qconsopt is public consumpition in levels, it DOES NOT take into account the degree of rivalry
        qconsopt = qconsopt0 
        zzpuopt = zzpuopt0
        utilopt = utilopt0  
        emaxopt = emaxopt0  
     ELSE
        stealopt = stealopt1
        newsavopt = newsavopt1
        vemaxopt = vemaxopt1  
        pconsopt = pconsopt1  
        !qconsopt is public consumpition in levels, it DOES NOT take into account the degree of rivalry
        qconsopt = qconsopt1 
        zzpuopt = zzpuopt1
        utilopt = utilopt1  
        emaxopt = emaxopt1  
     END IF
     
  END IF

!!$  IF(flag_here == 890) WRITE(*,'(a,8F20.10)') 'sol0', vemaxopt0, stealopt0, newsavopt0, pconsopt0, qconsopt0, zzpuopt0, utilopt0, emaxopt0
!!$  IF(flag_here == 890) WRITE(*,'(a,8F20.10)') 'sol1', vemaxopt1, stealopt1, newsavopt1, pconsopt1, qconsopt1, zzpuopt1, utilopt1, emaxopt1
!!$  IF(flag_here == 890) WRITE(*,'(a,8i5,10F20.10)') 'inps', time,imun,iedu,ipau,iwealth,ifu,ipr,iterm,ability,mwagein,wealthin,fundsin,zzprin,totrest,waget,lsmayor,R,netwealth
!!$  IF(flag_here==246) WRITE(*,'(a,9F20.10)') 'sol0 ', vemaxopt0, stealopt0, zzpuopt0, fundst, newsavopt0, pconsopt0, qconsopt0, utilopt0, emaxopt0
!!$  IF(flag_here==246) WRITE(*,'(a,9F20.10)') 'sol1 ', vemaxopt1, stealopt1, zzpuopt1, fundst, newsavopt1, pconsopt1, qconsopt1, utilopt1, emaxopt1
!!$  IF(flag_here==246) WRITE(*,'(a,9F20.10)') 'sol1 ', vemaxopt1, stealopt1, zzpuopt1, fundst, newsavopt1, pconsopt1, qconsopt1, utilopt1, emaxopt1

!!$  IF(flag_here==178) WRITE(*,'(a,1F50.10,8F20.10)') 'solsa', vemaxopt_sa, stealopt_sa, zzpuopt_sa, fundst, newsavopt_sa, pconsopt_sa, qconsopt_sa, utilopt_sa, emaxopt_sa

!!$  if (myrank == 0 .AND. flag_simulation == 1) write(*,'(a,10f20.10)') 'stealopt, pcons, qcons_w_riv', stealopt,pconsopt,qconsopt,qconsopt/((popsizeg)**eta),LOG(qconsopt/((popsizeg)**eta)),utilopt,emaxopt,vemaxopt,vemaxopt0,vemaxopt1

END SUBROUTINE mayordecision

FUNCTION focmayorTend(steal)
  USE commonvars
  IMPLICIT NONE
  REAL(8) :: focmayorTend, steal, pcons, zzpu, qcons, dqdpu, Esteal, fine_fun, Efine, finein
  INTEGER :: ifine

  EXTERNAL :: prodfunc, dprodfunc

  ! We compute expected stealing (taking into account the probability of being audited) instead of actual stealing as a trick to allow the differences in probability of being audited to have an effect at T
  
  ! Expected fine
  Efine = 0d0
  DO ifine = 1, Nfines
     frac_fine_g = frac_fine(ifine)
     finein = fine_fun(steal)
     Efine = Efine + finein*pfine(ifine)
  END DO
  Esteal = paudit_g*steal + (1d0-paudit_g)*(steal-Efine)

  IF (totrest + Esteal > 0d0) THEN
     pcons = totrest + Esteal
  ELSE
     pcons = small_no
  END IF
  !We use per-capita public inputs
  zzpu = (fundst - steal)/popsizeg
  IF (zzpu <= small_no) THEN
     qcons = subqcons
     dqdpu = 1000d0
     GOTO 670
  END IF

  !The output qcons is per-capita; we account for the degree of rivalry in the utility function
  CALL prodfunc(imunt,iedut,zzpu,zzprt,qcons)
  CALL dprodfunc(imunt,iedut,zzpu,zzprt,dqdpu)

  !We transform per-capita public consumption in total public consumption
670 qcons = qcons*popsizeg
!!$  dqdpu = dqdpu*popsizeg

  !We then transform the level of public consumption to take into account the degree of rivalry
  qcons = qcons/(popsizeg**eta)

  focmayorTend = pcons**(-gamma)-theta*(qcons)**(-delta)*dqdpu/(popsizeg**eta)
  
  IF (pcons <= 0.0d0 .OR. qcons <= 0.0d0) WRITE(*,*) 'CONSUMPTION IS NEGATIVE focmayorTend', pcons, qcons, steal, imunt,iedut,zzpu,zzprt

!!$  if (myrank == 0) write(*,'(A,9F30.10)') 'focmayor', focmayorTend, pcons**(-gamma), theta*(qcons)**(-delta)*dqdpu/(popsizeg**eta), pcons, totrest, Esteal, qcons, fundst, dqdpu
  
END FUNCTION focmayorTend

